#include <omp.h>
#include <iostream>
#include <math.h>

using std::cout;

double calcIntegral()
{
  
  const double minX = 0.0;
  const double maxX = 0.5;
  const double minY = 0.0;
  const double maxY = 1.0;
  
  const double step = 1e-3;
  
  int nx = static_cast< int >( ( maxX-minX ) / step );
  int ny = static_cast< int >( ( maxY-minY ) / step );
  
  double *xValues, *yValues;
  
  xValues = new double[ nx ];
  yValues = new double[ ny ];
  
  for( int i=0; i<nx; i++ )
  {
    double x = minX + i*step +step/2;
    xValues[ i ] = sin( x*M_PI ); 
  }
  
  for( int i=0; i<ny/2; i++ )
  {
    double y = minY + i*step +step/2;
    yValues[ i ] = cos( y*M_PI ); 
  }
  
  double sum = 0.0;
  
#pragma omp parallel for reduction(+:sum)
  for( int i=0; i<nx; i++ )
    for( int j=0; j<ny/2; j++ )
    {
      sum += exp( xValues[ i ] * yValues[ j ] );
      sum += exp( -xValues[ i ] * yValues[ j ] );
    }
    
  sum *= step * step;
  sum /= (maxX-minX) * ( maxY-minY );
  sum += 1;
  
  return sum;
  
}

int main( int argc, char **argv )
{
  
  omp_set_num_threads( 4 );
  
  double startTime = omp_get_wtime();
  cout << "Value: " << calcIntegral() << std::endl;
  double endTime = omp_get_wtime();
  
  cout << "Time: " << endTime-startTime << std::endl;
  
  return 0;
  
}